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We calculate radial migration rates of protoplanets in laminar minimum mass solar neb- 
^ ula discs using three-dimensional self-gravitating radiation hydrodynamical (RHD) models. 

The protoplanets are free to migrate, whereupon their migration rates are measured. For low 
mass protoplanets (10-50 M®) we find increases in the migration timescales of up to an or- 
der of magnitude between locally-isothermal and RHD models. In the high-mass regime the 
migration rates are changed very little. 
^ / These results are arrived at by calculating migration rates in locally-isothermal models, 

before sequentially introducing self-gravity, and radiative transfer, allowing us to isolate the 
eff'ects of the additional physics. We find that using a locally-isothermal equation of state, 
4!h without self-gravity, we reproduce the migration rates obtained by previous analytic and nu- 

Qh merical models. We explore the impact of diff'erent protoplanet models, and changes to their 

Q assumed radii, upon migration. The introduction of self-gravity gives a slight reduction of 

^ the migration rates, whilst the inertial mass problem, which has been proposed for high mass 

^ protoplanets with circumplanetary discs, is reproduced. Upon introducing radiative transfer to 

^ models of low mass protoplanets (^^ 10 M©), modelled as small radius accreting point masses, 

we find outward migration with a rate of approximately twice the analytic inward rate. How- 
^ ever, when modelling such a protoplanet in a more realistic manner, with a surface which 

^ enables the formation of a deep envelope, this outward migration is not seen. 
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1 INTRODUCTION mentum exchange was such that a planet's orbital radius should 
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^, ^ . ^ . . , , , . , change significantly on short timescales (~ 10^ - 10 years). 

The discovery oi giant planets m tight orbits around their stars has 

introduced a complication that planet formation models must ad- There are two primary types of planet migration which act in 
dress. The formation theories contemporary to the discovery of 51 diff'erent mass regimes. Type I acts in the low mass regime (< 100 
^ Peg b ( [Mayor & Queloz|1995| ), the first such exoplanet discovered, M®), and is a result of a low mass planet's interactions with an- 
0^ had assumed planet systems with morphologies similar to the solar nuli of gas that are in resonance with it, at so called Lindblad res- 
system, and were unable to explain the planet's creation. The high onances, as well as the corotation resonance. Angular momentum 
temperatures in the inner regions of circumstellar discs, as well as exchange occurs efficientiy with Lindblad resonances both inside 
the dominance of the star's gravity, prevent formation by gravita- and outside of the protoplanets orbital radius. The inner Lindblad 
tional instability. Fairing no better, the core accretion model is con- resonance passes angular momentum to the protoplanet, acting to 
founded by the low solids surface density in such high temperature accelerate it, whilst the outer Lindblad resonance works in oppos- 
regions where no ice can form. However, a mechanism that might tion. A linear analysis can be used to describe the eff'ect of Type I 
off'er an explanation for the existence of such Hot Jupiters had al- migration, and in doing so it is found that the Lindblad resonances 
ready been suggested by |Goldreich & Tremaine| ( |198 0') who had are not spaced symmetrically inside and outside of the planet's 
studied the rate of angular momentum and energy transfer between orbit. Instead each outer Lindblad resonance is somewhat closer 
a disc and an embedded planet. They found that the angular mo- than the corresponding inner resonance, allowing the angular mo- 



mentum exchange with these outer resonances to dominate ( Ward| 
|1986| ). The interaction of a planet with gas in the corotation region, 
the so called corotation torque, has a magnitude at most of only a 
* E-mail: ayliffe@astro.ex.ac.uk f^w tens of percent of the diff'erential Lindblad torque ( |Ward|1997{ 

t E-mail: mbate@astro.ex.ac.uk [Tanaka et al.|2002| ) when fully unsaturated. As a result, the net loss 
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of angular momentum by the planet to the disc due to the differ- 
ential Lindblad torque causes the planet to migrate inwards ( |Ward| 
[T986) . 

The rates of planet migration derived from analytic models 
(i.e. |Ward||1997[ [Tanaka et al.|2002] ) are in agreement with those 



obtained using locally-isothermal numerical models ( Korycansky 
& Pollack] 1993| [Nelson et al.|2000|passet,2002,.D'Angelo et al. 



2002[|Bate et al |2003[ |D'Angelo et al.|2003[|Alibert et al.|2004 



D^ngelo et al."'2005'; "Klahr & Kley"2006! 'D' Angelo & Lubow| 
2008 , Paardekooper & Mellema 2008a, Li et al. 2009 ). These rates 
are such that low mass planets very rapidly fall into their stars, on 
much shorter timescales than are required for them to grow to of or- 
der a Jupiter mass. As such, to explain the observed exoplanet pop- 
ulation, some modification is required to slow down Type I migra- 
tion, something beyond what the locally-isothermal models have 
considered ( |Alibert et al.|2004| [2005| [Ida & Lin|2008| |Mordasini| 
let al.|2009| ). 

One such modification is the inclusion of turbulence, dispens- 
ing with the usual laminar disc assumption. Three-dimensional 
magnetohydrodynamic models performed by |Nelson & Papaloizou| 
( |2004| ) focused on the implications for migration of interactions of 
protoplanets within magnetorotationally unstable accretion discs. It 
was found that the disc turbulence led to large scale fluctuations in 
the net torque experienced by a protoplanet over the course of its 
orbit. As such the change in a protoplanet' s orbital radius becomes 
a random-walk. The most significant fluctuations were found for 
protoplanets of < 30 M®, for which the time averaged fluctuations 
were inconsistent with conventional Type I migration. The increase 
in a low mass protoplanet' s migration timescale due to the oscil- 
lation of its orbit suggests turbulence as a route by which to im- 
prove the survival probability of such bodies dNelson & Papaloizou| 
|2004[ |Nelson|2005| ). Another addition to migration models which 
can lead to slower rates is more realistic thermodynamics. Sev- 
eral numerical models have introduced non-isothermal equations of 
state ('Morohoshi & Tanaka 2003 , Jang-Condell & Sasselov 2005 , 
[Paardekooper & Mellema 2006 , 2008b , Baruteau & Masset 2008a, 
Paardekooper & Papaloizou 2008b , Kley & Crida 2008 , Kl ey^etaf] 
|2009j), and have found that the rate of Type I migration can be re- 
duced, or its direction reversed. 

The second type of migration, called Type II, becomes eff'ec- 
tive once a planet grows to a mass of a few tenths of a Jupiter mass 
and is able to clear its corotation region of material, forming a gap 
in its parent disc (Ward 1997, Bry den etal.r i999, Kley 1999). The 
width of the gap is determined by the competing influences of the 
viscous and pressure forces, which tend to close the gap, and the 
gravitational torques which work to open it. The planet becomes 
locked in the gap it has created, and migrates inwards on the disc's 
viscous accretion timescale. Such Type II migration continues on 
the viscous timescale whilst the protoplanet mass is less than or 
comparable to the mass of the local disc with which it is interacting. 
A protoplanet that comes to dominate the local mass can no longer 
be treated simply as a constituent part of the disc. Rather, its inertia 
acts to slow the rate of migration that the disc's viscous evolution 
can bring about jSyer & Clarke|1995l|Ivanov et al.|1999| ). |Ivanov| 
|et al. I suggested an expression to describe this late stage migration 
which led to a reduction of the migration rate with increasing pro- 
toplanet mass and with reduced orbital radius. This gives a slowing 
migration that might enable massive protoplanets to avoid falling 
into the central star, and help to form Hot Jupiter s. The transition 
from Type I into Type II migration cannot be treated analytically 
due to the many non-linear processes at work, but is has been mod- 



elled numerically in both two and three -dimensions jBryden et"aL] 
|T999 Kley 19991 P'Angelo et al.|2002| [Bate et al.|2003l L 

This paper discusses numerical models conducted to investi- 
gate planet migration, as has been done previously by numerous 
authors in both two and three dimensi ons ([Korycan sky & Pollack 



1993^ 



N elson et al.'2000','Massetl2002',"D'Angelo et al. 2002, Bate 
etaLl2 003 , Nelson & Papaloizou 2003 , D' Angelo et al. 2003 , Alib- 
ert et al. 2004 , Nelson & Papaloizou 2004 , Nelson 2005 , D' Angelo 



et al, 



bidel 



2005, Mas set et al. 2006, Klahr & Kley 2006, Crida & Mor- 
|2007[ [Paa rdekooper & Papaloizou 2008b, Paardekooper & 
"vTTinnnoT r^'A^r^^i^ s-r t hU^^t, onnc i^i^^. s-r r^*.;^o onnoT 



Mellema'^2008a, D' Angelo & Lubow 2008, Kley & Crida 2008^ 
Peplinski et al. 2008a b , Li et al. 2009 , Paardekooper & Papaloizou 
2009, Crida et al. 2009, K ley etal.|2 009, Yu et al.|2010| [Baruteau 



,& Lin. 2010 Some of the more recent of these works have begun to 
include more physics, such as magnetic fields, and non-isothermal 
equations of state as discussed above. Other models, conducted by 
[Nelson et al. (200 0), D'Angelo erar|([2003]),[A libert et al.[([ 2004l ), 
and |Alibert et al.[p005| ), have included mass and angular momen- 
tum accretion onto migrating protoplanets, whilst the influence of 
self-gravity was explored by [Nelson & Benz[ ([2003]), [Baruteau &| 
[Masset[ ( [2508bt , and analytically by [Pierens & Hure[ ([2005). Re- 
gardless of the physics included, it has been common to exclude 
torques from the immediate vicinity of the planet due to the diffi- 
culty of resolving this region in any great detail. This poor resolu- 
tion also prevents the self-consistent modelling of a protoplanetary 
envelope, which requires both good resolution and the inclusion of 
self-gravity and radiative transfer. 

In this work we conduct three dimensional global disc mod- 
els, using smoothed particle hydrodynamics (SPH), that include ra- 
diative transfer, self-gravity, and which use a planetary surface to 
allow modelling of gas flow to well within the Hill sphere. With 
our surface treatment, and the adaptable resolution of SPH, the in- 
nermost torques can be well resolved and a realistic protoplanetary 
envelope can develop. We conduct models of varying complexity, 
introducing more physics in stages such that we can determine the 
impact of each addition on protoplanet migration. The first models 
we discuss are designed to emulate the treatment used in the grid 
models of [Bate et al.[ ( [2003] ), allowing us to test our SPH calcula- 
tions. We then introduce mass and angular momentum accretion, 
followed by the addition of self-gravity. Next the surface treatment 
is employed to explore the impact of realistic gas flow within the 
Hill sphere. Finally radiative transfer is folded into the mix, and a 
range of opacities explored. 

In Section 2, we describe our computational method and its 
testing. In Section 3 we present the results of our models, starting 
with the simplest locally-isothermal case, and progressing to self- 
gravitating radiation hydrodynamic cases using our planet surface 
treatment. Section 4 discusses the implications of our results for gi- 
ant planet formation theory, whilst a summary and our conclusions 
are given in Section 5. 



2 COMPUTATIONAL METHOD 

The calculations described herein have been performed using a 
three-dimensional SPH code. This SPH code has its origins in a ver- 
sion first developed by Benz[( [1990[[Benz et al.f l990 ) but it has un- 
dergone substantial modification in subsequent years. Energy and 
entropy are conserved to timestepping accuracy by use of the vari- 
able smoothing length formalism of [Springel & Hernquist[ ( [2002] ) 
and [Monagh an (2002J with our specific implementation being de- 
scribed in [Rice & Bate] ( [2007] ). Gravitational forces are calculated 
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and neighbouring particles are found using a binary tree. Radiative 
transfer is modelled in the two temperature (gas, Tg, and radiation, 
Tj-) flux-limited diffusion approximation using the method devel- 
oped by Whitehouse et al.|p005 ) and Whitehouse & Bate ( 2006] ). 
Integration of the SPH equations is achieved using a second-order 
Runge-Kutta-Fehlberg integrator with particles having individual 
timesteps ( |Bate|1995] ). The code has been parallelised by M. Bate 
using OpenMR 



2.1 Disc setup 

We model a protoplanetary disc (In radians) with radial bounds of 
0.1 - 3 Tp (0.52 - 15.6 AU), where is the initial orbital radius of 
our embedded protoplanets, taking a value of 5.2 AU. The radial 
temperature profile for the disc is initialised as Tg oc r~\ whilst 
the surface density of the disc goes as E oc r~^^^; these profiles 
are equivalent to those of [Bate et al.| ( [2003"] ), giving a constant ra- 
tio of disc scaleheight to radius of H/R ^ 0.05. At the initial 
temperature is ^ 75K and the surface density is 75 g cm~^. This 
leads to a total disc mass of ^ 0.005 M© within the boundaries, (or 
^ 0.0075 Mo within 1.56 - 20.8 AU inline with Bate et al. 2003 or 
0.01 Mo within 26 AU inline with D'Angelo et al. 2003 ). The disc 
begins with a Keplerian velocity structure, established around a 1 
Mo star. It should be noted that when modelled with radiation hy- 
drodynamics the disc settles to have a slightly steeper temperature 
profile at the midplane, where in our disc it be comes approximat ely 
Tg oc r~^^. A similar steepening is found by Kley et al. ( |2009| ) in 
their radiation hydrodynamical models. 

The inner edge of the disc is bounded by a potential barrier 
which acts to prevent the disc accreting on to the central star. The 
outer edge of the disc is surrounded by a population of ghost parti- 
cles that act as a pressure barrier to prevent the disc from spreading 
due to shear viscosity. Particles that do move beyond the outer edge 
of the computational domain are removed from the calculation. The 
migration rates of planets are determined by the structure of the disc 
in their vicinity, particularly the surface density profile. As such, 
the aim of these boundaries was to establish a broad region of disc 
about Tp throughout which the surface density followed the desired 
profile, with the only modification being due to the planet-disc in- 
teraction. Maintaining a broad stable region is important to allow 
any embedded protoplanet to migrate away from its initial position 
at Tp, and to allow density waves to propagate to large distances, 
such that their interactions with the protoplanet are included to the 
point of inconsequence. 

The inner boundary allowed gas to pile up at a small radius, 
far from the intended initial orbital radius of any introduced pro- 
toplanets. Once gas has accumulated at the boundary, the concen- 
tration of material exerts an outward pressure gradient back into 
the disc, supporting further gas against inward migration. Fig.[T]il- 
lustrates the disc profile that develops, and illustrates the success 
of the boundaries in preserving the desired surface density profile 
(shown by the long-dashed line) throughout the region of likely 
migration, from 0.5 to 2 (2.6 - 10.4 AU). The profiles shown 
are for calculations with 5 x 10^, 2 x 10^ and 10^ particles. In the 
highest resolution calculation, the rarefaction due to the outwards 
pressure from the boundary is narrowest and shallowest. However, 
the radial range over which a consistent surface density profile is 
achieved is not greatly changed between the two higher resolution 
models. For calculations performed with 5x10^ particles this is not 
the case, with a much greater encroachment of the boundary eff'ects 
into the likely migration region. In the calculations that follow we 
use 2x10^ particles as this was sufficient to establish a stable disc 



100 




Figure 1. Surface density profiles of a locally-isothermal disc modelled 
with lO'^ (solid line), 2x10^ (short-dashed line), and 5 x 10^ (long-dashed 
line) particles, in the absence of a planet, after 4 orbits of the disc's outer- 
edge. The straight dashed line across the top denotes the desired surface 

_ 1 _^ 
density profile of E(r) oc r 2 , with a value of 75 g cm at r = rp . With the 

implemented disc boundaries this profile is achieved by the two higher res- 
olution calculations in the region over which a planet may migrate, whilst 
in the lowest resolution calculation the eff'ect of the inner boundary reaches 
to almost Tn. 



in the planet's vicinity, whilst being considerably faster than using 
10^ particles, for which the benefits were marginal. 

Upon creating a disc model, there is generally a period of tran- 
sient structural perturbations as the model settles to a stable state. 
This settling must be completed before protoplanets are introduced 
in order that their orbits are not immediately disturbed by perturba- 
tions unrelated to the interesting protoplanet-disc interactions. The 
initial discs were evolved in the absence of a planet until any tran- 
sience resulting from settling had dissipated, which required just 
over 4 orbits of the disc's outer edge. 



2.2 Equation of state 

Two equations of state are used for the calculations presented in 
this paper. One is a locally-isothermal equation of state with the 
temperature of the gas throughout the disc being a fixed function of 
radius (as described in the disc setup). The second is used for our 
radiation hydrodynamical calculations, which are conducted using 
an ideal gas equation of state p = pTgRg//u where Rg is the gas con- 
stant, p is the density, Tg is the gas temperature, and is the mean 
molecular mass. The equation of state takes into account the trans- 
lational, rotational, and vibrational degrees of freedom of molecu- 
lar hydrogen (assuming a 3:1 mix of ortho- and para-hydrogen; see 
[Boley et al.|2007 l). It also includes the dissociation of molecular hy- 
drogen, and the ionisations of hydrogen and helium. The hydrogen 
and helium mass fractions are X = 0.70 and Y = 0.28, respec- 
tively, whilst the contribution of metals to the equation of state is 
neglected. More details on the implementation of the equation of 
state can be found in [Whitehouse & Bate| ( |2006| ). 
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The two temperature (gas and radiation) radiative transfer in 
the flux-Umited diffusion approximation employed in this work is 
implemented as described by jWhitehouse et aL] ( |2005| ) and [White] 
[house & Bate ( 2006) . Briefly, work and artificial viscosity (which 
includes both bulk and shear components) increase the thermal en- 
ergy of the gas, and work done on the radiation field increases the 
radiative energy which can be transported via flux-limited diff'u- 
sion. The energy transfer between the gas and radiation fields is 
dependent upon their relative temperatures, the gas density, and the 
opacity, k. 



2.3 Opacity treatment 



We use interpolation of the opacity tables of jPollack et al.| ( p^85] ) to 
provide the interstellar grain opacities for solar metallicity molec- 
ular gas, whilst at higher temperatures when the grains have subli- 
mated we use the tables of |Alexander| ( |1975] ) (the IVa King model) 
to provide the gas opacities (for further details see | Whitehouse| 
|& Bate,2 006). The grain opacities are based on interstellar abun- 
dances, and may not hold true in a circumstellar disc where ag- 
glomeration can strongly modify the grain population, reducing 
the opacity. Therefore, we vary the grain opacity in our models to 
mimic the eff'ects of these processes; the gas opacities of | Alexander] 
( p^5 ) that are used at higher temperatures are not scaled. In this 
paper we use 100, 10, and 1 per cent interstellar grain opacities (see 
[Ayliff^e & Bate|2009a| for further details). 



2.4 Radiation boundary 

Many of the calculations performed here involved radiative trans- 
fer, and so required a mechanism by which energy could be lost 
from the disc into the encompassing vacuum. This necessity arises 
because the flux-limited difl'usion scheme transfers energy between 
SPH particles, rendering it unable to radiate into a vacuum where 
there are no particles. By determining the height above and below 
the midplane at which the optical depth into the disc. Top ~ 1, we 
can identify two layers of particles which exist in optically thin re- 
gions. Compelling particles that fall within these regions to follow 
the initial temperature profile of the disc allows the energy they re- 
ceive from the bulk of the disc to be eff'ectively lost, as though into 
the vacuum. This initial temperature profile was chosen for pur- 
poses of comparison with previous work where it is common, and 
crudely represents the insolation of the disc. 

The vertical density structure of the disc approximately fol- 
lows a Gaussian distribution, integration over which yields the sur- 
face density, E^, in terms of the error function, erf(). The optical 
depth is given by Top (a) = - ^ xpdz. Completing this integration 
with the Gaussian form of p, setting Top = 1, and the upper limit of 
integration set to infinity to represent the vacuum, it is possible to 
rearrange to find the lower limit of integration, depth 'a', as shown 
in equation [T] 



a 
H 



V2erf-Ml 



2top 

KLr 



(1) 



This is the depth in the disc from above which radiation can escape 
to the vacuum, and so defines the height of our radiation boundary. 
The boundary height varies with radius, due to the radial depen- 
dence of and /c, the values of which are taken from the initial 
disc distribution. 



2.5 Planet treatments 

Once a suitable disc has been established, it is possible to embed a 
protoplanet and model its evolution. In all of the models discussed 
in this paper, the protoplanet is free to move from its initial orbital 
radius, and its migration rate is measured, rather than calculated 
from static torques. How the protoplanet interacts with the disc, 
and so migrates, depends upon the way in which the protoplanet 
is modelled, and we employ several methods in this work. In the 
first case gas is not accreted, but instead gas that comes within 
a specified radius, race, of the point mass representing the proto- 
planet, is removed from the calculation, and its energy, momentum, 
and mass are all dispensed with. We call the particles used for this 
type of protoplanet treatment Killing sink particles, and they mi- 
grate purely due to torque interactions with the disc, approximately 
recreating the treatment used in [Bate et"ar] ( |2003j . Typically we set 
race to be some fraction of the Hill radius, R^. 

The second treatment also models the protoplanet as a sink 
particle which removes gas that comes within race, but in this case 
the gas properties are added to the sink particle. These sink parti- 
cles are called Accreting sinks. They are able to absorb the energy, 
mass, linear momentum, and angular momentum of the accreted 
gas (see |Bate et al.|1995| ). 

The third and final protoplanet treatment is one that we have 
used in previous work ( [Aylifi'e &^a te 2009a b ) which uses a sur- 
face force that allows gas to pile up on the surface of a protoplanet, 
realistically modelling accretion. This scheme is now applied by 
wrapping the surface around a gravitating point mass that is free 
to move, and which represents the protoplanet's initial mass. Such 
a protoplanet can move in response to changes in the momentum 
of the point mass and its bound envelope, the extent of which is 
self-consistently arrived at, rather than being prescribed. The sur- 
face is applied as a modification to the usual gravitational force in 
the protoplanet's immediate vicinity, taking the form 



Fr 



1 - 



2Rn 



(2) 



for r < 2 /?p, where is the protoplanet radius, r is the radius 
from the centre of the protoplanet, and Mp is the protoplanet's ini- 
tial mass. This equation yields zero net force between a gas parti- 
cle and the gravitating point mass at the surface radius R^. Inside 
of this radius the force is outwards and increases rapidly with de- 
creasing radius, such that the inner most gas particles come to rest 
very close to the surface, with an equilibrium position slightly in- 
side due to the envelope weight pressing down. As in our previous 
work, we provide a smooth start to the calculations by embedding a 
protoplanet of radius R^ - 0.0 Irp which then shrinks exponentially 
to the desired radius during the first orbit of the protoplanet. 



2.6 Resolution dependence of migration 

Resolution was considered previously in relation to the disc bound- 
aries, and establishing the desired density profile in an undisturbed 
disc. However, ensuring that the planet-disc interaction is well re- 
solved in our models, and that the migration rates obtained are in- 
dependent of resolution, required further testing. To this end, two 
simulations were performed which were identical in every respect 
excepting their resolution. A 10 M© protoplanet, modelled with an 
Accreting sink particle, was used for these calculations. The spiral 
waves launched by such a low mass planet are weak, as are the re- 
sulting torques acting back upon the planet. So whilst a degree of 
migration is expected, insufficient resolution might poorly resolve 



© 0000 RAS, MNRAS 000, 000-000 



Planetary migration and growth 5 




^ 0.996 



6 



20 30 
Time [Orbits] 

Figure 2. The orbital evolution of a 10 M© protoplanet modelled using 
Accreting sink particles, with accretion radii of -0.2 Ru (a radius chosen 
to mimic the ZEUS models of B ate et al.|20 03 ). Two resolutions were used 
to perform otherwise identical calculations, 10^ particles (thick upper line) 
and 2x10^ particles (thin lower line). The evolution covers 50 orbits of the 
protoplanet using a locally-isothermal equation of state. The model at both 
resolutions establishes the same rate of migration after a period of settling, 
roughly equal to the libration time. 



these small effects; this makes such a planet an ideal candidate for 
testing resolution. Fig.[2]shows the evolution of the planet's orbital 
radius with time for two resolutions. 

There is an initial period of settling whilst the torques due to 
the horseshoe region are established. As was seen in |Masset|p002| ), 
the outward corotation torques almost balance the inward differen- 
tial Lindblad torques when the protoplanet is first embedded lead- 
ing to almost no migration. However, as the model evolves the 
corotation torques become partially saturated, reducing their mag- 
nitude and allowing the Lindblad torques to dominate. As a result 
the protoplanet' s migration is seen to accelerate until the corota- 
tion torques reach a quasi-equilibrium state, which takes of order a 
libration time ( |Kley et aL|2009j ), given by 



2.7 Viscosity 

Viscosity is the means by which angular momentum is transported 
through a laminar circumstellar disc, enabling material to flow in 
towards the star. The viscosity envisaged is not caused by a molec- 
ular interaction, which is insufficient to explain the rates of accre- 
tion inferred from observations. The primary alternative is thought 
to be a magneto-rotational instability ( ,Balbus & Hawley|1991| ). 

Within the models presented here the viscosity is not achieved 
by modelling the responsible processes, but instead by using a pa- 
rameterised form of viscosity, though somewhat different to a typ- 
ical Shakura-Sunyaev of- viscosity ( [Shakura & Sunyaev|1973| ). An 
artificial viscosity for SPH, designed to conserve linear and angu- 
lar momentum, was introduced by |Monaghan & Gingold| ( |1983| ), 
and modified by Monaghan (1992 ) to deal with high mach num- 
ber shocks. This artificial viscosity is not designed to represent real 
viscosity, but instead to allow the modelling of shock phenomena, 
and to damp numerical noise. The a and terms set the strength 
of the viscosity, with typical values of a = 1 and p - 2. The two 
terms serve different functions. The or-term establishes a viscosity 
to damp subsonic velocity oscillations that may be produced in the 
wake of a shock front. The /?-term prevents particle interpenetration 
in supersonic shocks. 

The viscosity is only applied when the gas is under compres- 
sion, ideally near shocks, and we use the switch developed by |Mor-| 
|ris & Monaghanj ( p^ 7 ) to try and reduce the action of artificial 
viscosity where the cause is not a shock. Using this switch means 
that instead of the global a being applied to all the particles, they 
each scale this value (down to a minimum value of ofmin) based on 
their proximity to a shock. Setting = 0.1 can significantly re- 
duce unwanted dissipation. The differential rate of rotation of the 
disc leads to a shear viscosity by which angular momentum is ex- 
changed through the disc. The value of the viscosity at any point 
can be stated as a Shakura-Sunyaev a-viscosity, c^ss, using 



H 



(4) 



as given by |Lodato & Price| ( |201 0), where the coefficient given here 
is a factor 2 smaller because the viscous force is only applied be- 
tween approaching, not receding, particles. Here h is the local aver- 
age SPH smoothing length, and H is the disc scaleheight. Around 
Tp these quantities have values of approximately 0.016, and 0.05 
respectively, with a ^ 0.25 giving ofss ~ 0.004. This is typical for 
such a disc, and inline with the models of , Bate et al.) ( |2003j upon 
which our initial conditions are based. 



T//Z,= ^-Pp, (3) 
3 Xs 

where Xs is the width of the horseshoe region, and Pp is the planet's 
orbital period (11.8 years in the calculations presented here). Fol- 
lowing Kley et al. by using an approximation for the width of 
the horseshoe region, taken from the locally-isothermal models of 
[Masset et aL| ( |2006| , it is possible to get an idea of the horseshoe 
settling time. For a 10 M© protoplanet Tm, ~ 50 orbits, whilst for 
a 333 M© protoplanet Tm, ~ 8. It appears from Fig. |2] that a clear 
migration trend is established in around 30 orbits. Once Type I mi- 
gration is established, the migration rate is consistent between the 
two different resolution calculations. It is not clear why the initial 
settling appears to be resolution dependent, but the agreement of 
the eventual rates gives us confidence in the migration rates ob- 
tained using 2x10^ particles. 



2.8 Measuring the migration rate 

The migration rates were calculated by taking linear fits of the or- 
bital radii evolution data (e.g. Fig. [2]). Each simulation was allowed 
to evolve for at least 50 orbital periods, (as was done in [Bate et aL| 
|2003p , of which only the latter 25 were used for fitting purposes. 
The early orbits are sacrificed to ensure that any transient disrup- 
tion to the disc caused by the sudden introduction of a planet have 
subsided and to allow gas in the horseshoe region to reach a quasi- 
equilibrium state, as discussed previously. Fig. [3] illustrates the or- 
bital migration of 4 different mass protoplanets modelled with Ac- 
creting sinks and a locally-isothermal equation of state. It can be 
seen that the final 25 orbits of evolution are settled, and so provide 
a suitable period over which to measure the migration rates. For the 
highest mass protoplanets it can be seen that their migration over 
the last 25 orbits is linear in time, suggesting that the rate is well 
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Figure 3. The orbital evolution of 4 protoplanets modelled in locally- 
isothermal discs, using Accreting sink particles with race = 0.1 /?h- In all 
instances the planetary migration is well underway and free from large scale 
fluctuations after 25 orbits. The migration rates are measured over the last 
25 orbits. The Jupiter mass protoplanet appears to settle in ^ 10 orbits, 
inline with the libration time of 8 orbits jKley et al.|2009) . 



established. Deviations from linearity, which might cause a long 
term change in the rate are very subtle, and would only be estab- 
lished on very long timescales, inaccessible to our models due to 
their computational expense. 

There are instances where measuring the migration rate over 
the last 25 orbits will lead to faster rates for a stated mass than 
should truly be ascribed to it. In the first 25 orbits, in the Ac- 
creting sink and point mass with surface calculations, the proto- 
planet' s mass has increased through accretion. As a result of this, 
the migration rate tends to accelerate, giving a non-linear orbital 
radius evolution. It is therefore the case that linear fits to the latter 
half of this evolution don't exactly render the gradient expected for 
the initial protoplanetary mass. The impact of such mass accretion 
most strongly influences the orbital evolution of the lowest mass 
cores, which accrete a higher fraction of their total mass. As a re- 
sult, we present the migration rates with uncertainties in the mass, 
portrayed by error bars which stretch across the mass range over 
which the linear fit is taken, with the migration timescale plotted at 
the mean mass. In the Accreting sink calculations, the mass of the 
protoplanet includes that of all the gas accreted, but for the point 
mass with surface calculations the accreted mass is measured as 
the mass within the Hill radius. For a majority of the calculations, 
the mass spread is very small, and thus the mass associated with 
a given migration rate is well known. The migration timescale, r, 
is then calculated as the time a planet of fixed mass would take to 
migrate from its starting orbital radius, (5.2 AU in all of the models 
discussed), in to the central star; r = rp/r. The rate of migration 
is assumed constant in calculating these timescales, allowing ready 
comparison with previous works. 



10-^ 



Mp [Earth masses] 
10 100 



1000 



106 



105 



104 



Ward 1997 



/ \ 

/ 



\Tanaka et al. 2002 
VlJ 



0.01 0.1 1 

M [Jupiter masses] 



10 



Figure 4. A reproduction of figure 10 from [Bate et al.|p003) showing the 
migration rates they obtained by considering the torques at work upon pro- 
toplanets of fixed orbital radius in their models. The dots (with solid error 
bars) denote the rates calculated when the torques outside of Ru were in- 
cluded, the lower error bars are the rates when including torques external 
to 0.5 Ru, and the upper error bars are the rates for just the torques beyond 
1.5 Rn. The analytic Type I model of |Tanaka et"n^ ( p002) and the Type 
I/II model of Ward |T997j are plotted as short-dashed and long-dashed lines 
respectively. The triangles (with dotted error bars) mark the measured mi- 
gration timescales from the SPH models in which the protoplanet is able to 
migrate; the Killing sink particle accretion radii are Ru- The lower and up- 
per error bars for the SPH cases mark the measured rates with Killing sink 
particles of 0.5 and 1.5 Ru accretion radii respectively. The SPH models 
show good agreement with the grid-based ZEUS models. 



3 RESULTS 
3.1 Killing sinks 

We began by demonstrating that out SPH code can reproduce the 
migration rates obtained from past studies. These comparisons are 
between the migration rates obtained using our Killing sinks in 
locally-isothermal discs, and those obtained in ZEUS finite dif- 
ference models by Bate et al.| ( |2003| ), who used the same initial 
conditions. The Killing sinks are used to approximately imitate the 
planet treatment of 'B ate et al.| We expect our SPH models to re- 
produce similar migration timescales to these previous grid based 
results, however. Bat e et al. [ calculate their migration timescales by 
examining the net torque due to gas outside the planet's Hill radius 
upon the planet, which was on a fixed orbit at 5.2 AU. In making 
this comparison, the implicit assumption is that the torques do not 
change significantly when the planet is allowed to migrate rather 
than being held on a fixed orbit. Indeed, over the 50 orbits that are 
modelled in the SPH calculations, the planets do not migrate any 

A reproduction of 
4l where the ZEUS 



considerable distance, as can be seen in Fig. p 
Fig. 10 from jBate et al.| ( [2003] ) is shown in Fig. 
results are plotted with circles and solid error bars, whilst the SPH 
derived migration timescales are plotted using triangles with dotted 
error bars. The term error bar is a misnomer, instead the extremes 
of these bars indicate the migration timescales calculated by in- 
cluding more (0.5 Ru, lower bound) or less (1.5 Ru, upper bound) 
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Figure 5. A comparison of the protoplanet migration timescales measured in locally-isothermal, non-self-gravitating models using Killing sinks (asterisks), 
that exclude mass and angular momentum accretion, and otherwise identical calculations using Accreting sink particles (plus symbols) which include such 
accretion. The left panel is for those calculations where race = 0-5 Rn, whilst in the right panel race = 1.0 /?h- The analytic Type I model of Tanaka et al. (2002/ 
is plotted as a dashed line. Note that the Accreting sinks begin with the same four masses as the Killing sinks, but their masses increase substantially during 
the calculations, particularly the low-mass protoplanets. The mass error bars for the Accreting sink particles show the range in mass over which the linear fit 
to determine the migration rate was applied. Including mass and angular momentum accretion consistently leads to faster rates of migration when measured at 
equivalent masses. 



of the gas closest to the protoplanet in the torque calculations. In 
the SPH models these same limits are found by conducting simu- 
lations with killing radii of race = 0-5 Rn^ 1 Rn, and 1.5 /?h; note 
that this is not exactly equivalent to changing the torque exclusion 
radius, and the limitations of the comparison will be discussed in 
the following paragraphs. For all protoplanet masses the migration 
timescales obtained with the SPH code, using killing radii of /?h, 
differ by less than 35 per cent from the ZEUS based results. The 
SPH results also show the expected transition from the linear Type 
I migration regime to the non-linear Type II beyond ~ O.lMj. As a 
result, we are confident that our SPH code is capable of modelling 
migration proficiently. 

The spread in the migration timescales we find for the various 
killing radii are larger that those found in the ZEUS calculations 
with diflferent torque exclusion radii. Explanations for these larger 
spreads can be found in the density and resolutions of the two mod- 
els. For a given mass protoplanet, there is only one ZEUS model 
upon which the three diflferent torque exclusion radius calculations 
are based. This gives a consistent density structure between the dif- 
ferent cases, but a fixed spatial resolution which is poor near to 0.5 
/?H- In the SPH calculations, the race equivalent to these exclusion 
radii are each probed in separate models. The evacuated region for 
an race = 1.5 Killing sink is far larger than for the 0.5 R^ case, 
which changes the circumplanetary density structure between the 
two cases, with lower densities at equivalent radii for the larger race- 
This leads to lower torques acting from the boundary of a Killing 
sink with race = 1.5 /?h than at an equivalent exclusion radius in the 
ZEUS models; hence the larger migration timescale. 

For race = 0.5 /?h the SPH resolution is far better than that 



of the ZEUS calculations, and so the models better discern the gas 
that is available to exert torques, leading to faster migration. At all 
the radii, regardless of the diflferent spreads, there is consistently a 
turn away from Type I to Type II migration at higher masses, and 
the absolute migration timescales are within a factor of two or bet- 
ter at each mass and radius. The figure makes clear that in both the 
SPH and ZEUS models, the migration timescales reduce when the 
killing radius (or torque exclusion radius) is shrunk. This indicates 
that the contributions from torques from within /?h are negative, as 
was found in |Bate et al.| ( |2003| ). This is an important agreement as 
the region within /?h is not treated in analytic descriptions of mi- 
gration, and so is not well understood. It appears that the persistent 
torques from within R^ are negative, which is in the opposite sense 
to the corotation torques in a disc with the surface density profile 
modelled here. We return to this in Section [34l 



3.2 Accreting sinks 

Whilst Killing sinks provided a means by which to test and com- 
pare our models, they do not allow us to include protoplanet 
growth. To that end we next introduced accretion into the migration 
models, using our Accreting sinks. Fig. [5] compares the migration 
timescale found for Accreting sink particles in locally-isothermal 
models with those obtained using the Killing sink method pre- 
viously described. Whilst the Killing sinks maintain their initial 
masses throughout, the Accreting sink models lead to increasingly 
massive protoplanets, sometimes growing to several times their ini- 
tial mass. This makes interpreting the relative migration rates some- 
what difficult, though the broad result is that including accretion 
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leads to faster migration rates. Comparing an Accreting sink which 
has grown from an initial mass of 10 M© to beyond 33 M© with a 
Killing sink which maintains a mass of 33 M© throughout its evo- 
lution must be done with some care as there are several factors 
contributing to their differing migration timescales. One factor is 
that the Accreting sink has a smaller r^^c than the 33 M© Killing 
sink model as the former is calculated for its initial mass of 10 M©, 
and so it develops a different circumplanetary structure. As was dis- 
cussed in section [3T| for the different sized Killing sinks, a smaller 
race leads to faster migration. This factor therefore contributes some 
fraction of the difference in migration rates seen in Fig. [5] 

As seen in Fig. |4] and again in comparing the left and right 
panels of Fig. |5] increasing race reduces the migration rates of 
the Killing sinks. However, the Accreting sink migration rates are 
largely unaffected by the change in race from 0.5 to 1 R^. The 
torques experienced by these sinks are equivalent to those of the 
Killing sinks. It therefore must be the case that the loss of nega- 
tive torques from within as race is increased is countered in the 
case of Accreting sinks by greater mass accretion. To test this we 
performed a series of calculations in which Accreting sinks were 
embedded in discs with which they did not gravitationally inter- 
act. These sinks accrete only material that they sweep up from their 
path through the disc as they cannot gravitationally draw in gas. 
This also means that they experience no torques from the disc. Any 
migration of these protoplanets is a result solely of the dilution of 
their specific angular momentum due to accretion of pressure sup- 
ported, and so sub-Keplerian, gas. 

To test this we numerically integrate from a protoplanet's ini- 
tial mass and radius to its final mass, accounting for its accretion 
of gas with sub-Keplerian angular momentum to compare its final 
orbital radius with that obtained from the SPH models. In our un- 
perturbed disc gas, at the initial orbital radius of the protoplanet, 
the specific angular momentum of the gas has a value of 0.998 of 
the protoplanet's specific angular momentum. However, we note 
from the SPH models that the mean radius from which gas is ac- 
creted is smaller than rp, and so the average specific angular mo- 
mentum is somewhat smaller than this value. Employing this mea- 
sured value in our simple integration we obtain final protoplanet 
orbital radii that are in good agreement with those measured from 
the no-interaction SPH models. Applying this iterative calculation 
to the Accreting sink models shown in Fig. |5] can account for a 
surprisingly large fraction of the change in their orbital radii. For 
a 10 M© protoplanet, about 50% of the change can be attributed 
to the effects of accreting sub-Keplerian material. This implies up 
to a factor of 2 difference in the migration timescale between an 
Accreting sink and a Killing sink due to gas accretion. 

Fig.[6]includes all the different accretion radii used in our Ac- 
creting sink calculations. It can be seen that those with radii of 
0.1/?H (marked as plus symbols) and masses < 0.3 Mjupiter closely 
follow the analytic model of |Tanaka et al.| ( [2002| ). In these cases 
the dilution of the angular momentum and the effective change in 
the ratio of race to /?Hare minimal due to the small accretion radii, 
and lower accretion rates. The accretion rates onto these Accreting 
sinks are still substantial because of the evacuation of gas within 
race, which leads to an inward pressure gradient near their bound- 
ary where the gas is removed. This effectively pulls material in at a 
faster rate than would be expected in the absence of the hole. This 
artificial evacuation at the planet boundary not only leads to artifi- 
cially fast accretion, it also alters the density structure around the 
protoplanet out to distances that are large compared to race; this is 
tackled in section [34l 
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Figure 6. Migration timescales for 10, 33, 100, and 333 M© protoplan- 
ets undergoing mass and angular momentum accretion. Each protoplanet is 
modelled as an Accreting sink particle, with race = 0.1 /?h (plus signs), 0.5 
Rw (diamonds), or 1.0 Rw (asterisks) for each mass. The analytic Type I 
model of |Tanaka et al.| ((2002 ) is plotted as a dashed line. The smallest ac- 
cretion radii protoplanet models give the migration rates that most closely 
match the analytic model. 



3.3 Self-gravity 

All of our models discussed to this point have included the gravi- 
tational interaction of the protoplanet with the disc and vice versa, 
but have neglected the disc's self-gravity, that is its influence upon 
itself. For a low mass disc this is a reasonable approximation as 
the central star is so dominant. Indeed the inclusion of self-gravity 
in low-mass discs has been found to have only a small impact on 

Pierens| 



migration rates by several authors ( Nelson & B enz 20031 
|& Hure|2005l|Baruteau & Masset|2008b[|Crida et al.|2069| . Com- 
paring two series of identical calculations, with the inclusion of 
self-gravity being the only differentiating characteristic, we found 
there to be a very small, but consistent increase in the migration 
timescales for self-gravity models. Fig.|7]makes evident the small 
magnitude of the change due to self-gravity, in this instance for Ac- 
creting sinks with race = 0.1 i?H- However, self-gravity remains an 
essential component in building a self-consistent model. 

When substantial circumplanetary discs or envelopes form, 
the inclusion of self-gravity becomes important in delivering ac- 
curate migration rates. As discussed in Crida et ah] ( |2009 ), a cir- 
cumplanetary disc is locked to the protoplanet, and moves with it 
around the central star. In the absence of the disc self-gravity the 
circumplanetary disc does not experience the torques due to the 
protoplanetary disc; the torques which the protoplanet is respond- 
ing to. As such it does not migrate of its own accord, but must be 
pulled by the protoplanet, artificially slowing the rate of migration; 
known as the 'inertial mass problem' ( |Crida et al.|2009| ). To explore 
the impact of this phenomenon we performed two calculations us- 
ing point masses with surfaces (discussed in the next section) to 
model 333 M© protoplanets, one model with and one without the 
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Figure 7. Comparisons between the migration timescales of protoplanets, 
modelled in locally-isothermal discs by Accreting sinks with r^cc = 0.1 
Ru, in non-self-gravitating discs (plus symbols), and self-gravitating discs 
(asterisks). As has been found by several authors, the inclusion of disc self- 
gravity has a very small effect upon the migration timescales. This small 
impact is consistent across all protoplanet masses, slightly increasing the 
migration timescales. The analytic Type I model of ,Tanaka et al.,^2QQ2^ is 
plotted as a dashed line. 



inclusion of disc self-gravity. The two resulting migration evolu- 
tions are shown in Fig. [8] In both cases a protoplanetary disc devel- 
ops, and in the case without disc self-gravity, this causes the pro- 
toplanet' s rate of migration to slow by approximately 14 per cent. 
Having illustrated the effect, we proceed with calculations which 
all include self-gravity to avoid what is an unphysical interaction. 



3.4 Protoplanets with surfaces 

The final addition to our locally-isothermal models was to use our 
planet surface treatment. This addresses the failing of our Accret- 
ing sink models by allowing the gas in the protoplanet' s vicinity to 
develop a structure free of the undesired effects of evacuating gas 
near the protoplanet. Combining this treatment with the inclusion 
of self-gravity we obtain a model in which the interaction of the 
disc and protoplanet is free of any artificially imposed radial limit, 
such as excluding material within R^. Instead, the radius at which 
gas becomes bound to the protoplanet, and stops applying orbit af- 
fecting torques is achieved self-consistently. Fig. [9] illustrates the 
midplane density distribution around both an Accreting sink and 
a point mass with surface (in both cases a 333 M© protoplanet, of 
radius 0.1 Ru). Whilst gas locked in a circumplanetary disc will 
not exert migration causing torques, gas beyond such a disc (found 
to be >Rn/3 in these models, [AylifFe & Bate 2009b') that moves in 
and out of the Hill sphere will. The density distribution in the two 
cases are not similar until radii of > /?h. suggesting that the Accret- 
ing sinks experience smaller torques from within and around Ru 
than the point masses with surfaces. The importance of these local 
torques can be seen in the left panel of Fig. [T0| which is a map of the 
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Figure 8. The orbital evolution of a 333 M© protoplanet modelled using a 
point mass with surface of radius 0.03 Ru, using a locally-isothermal equa- 
tion of state. A calculation including disc self-gravity is shown using a nar- 
row line, whilst the non-self-gravity calculation is shown using a thick line. 
There is reduction in the migration rate of the protoplanet modelled with- 
out self-gravity as it develops a circumplanetary disc; a result of the inertial 
mass problem. The migration rate is slower by 14 per cent. 
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Figure 9. Midplane densities around a 333 M© protoplanet modelled us- 
ing a locally-isothermal equation of state and self-gravity after 50 orbits. 
The solid line denotes the densities around a protoplanet modelled using a 
point mass with a surface of radius 0.1 Ru. The dashed line gives the mid- 
plane densities around an identical protoplanet modelled by an Accreting 
sink with race = 0.1 Ru- The evacuation of material at the Accreting sink's 
boundary leads to an inward pressure gradient that pushes in gas, giving 
artificially high accretion rates, and reducing the density in the vicinity of 
the growing planet out to radii of > /?h- 
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Figure 10. Torque distributions around a 333 M© protoplanet modelled by a point mass with a surface of radius 0.03 Ru after 50 orbits. The left hand panel 
is from a locally-isothermal model, whilst the right hand panel is from a radiation hydrodynamic calculation using interstellar opacities; both cases include 
self-gravity, and the figures exclude torques from within RuP- The Roche lobe is marked on in each case using a white line. All the torques above the white 
dashed line are positive, whilst all those below are negative. Introducing radiative transfer smears out the spiral arms, fills in the corotation region, and reduces 
the maximum torques which act from within the Roche lobe. The section shown is rp ± 8 Ru, which avoids the inner and outer boundary regions. 



torques acting on the protoplanet in the plane of the disc. The fig- 
ure excludes torques acting from within RuP which is the expected 
limit of the circumplanetary disc for such a high mass protoplanet, 
but note that within the model only the radius of the planet's surface 
is specified. Between RuP and the Roche lobe (outlined in white) 
are the strongest interactions by virtue of the gas's proximity and 
density. These torques are enabled by gas flow into and out of the 



Roche lobe which allows for the continuous exchange of angular 
momentum, establishing a persistent torque. This flow can be seen 
for locally-isothermal models in the top panels of Fig.[TT] The ve- 
locity vectors are plotted over a midplane density map of the area 
surrounding the protoplanet. Gas can be seen passing into and out 
of the Hill sphere on trajectories that are both radial (horseshoe or- 
bits) and azimuthal. 
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Figure 11. Sections of the whole disc models plotted as midplane density maps, focussing on the region surrounding the embedded protoplanets, with velocity 
vectors overplotted. It can be seen that material is flowing into the Hill sphere and back out again, reaching significant depths within this region. This material 
is able to exchange angular momentum with the protoplanet before escaping its potential, and so a persistent torque is established from within Rn. From left 
to right the protoplanet masses are 33, 100, and 333 M©. From top to bottom the models are locally-isothermal, and then radiation hydrodynamical with 1, 10, 
and 100 per cent interstellar grain opacities. 



Fig. [12] includes migration timescales from Accreting sink 
models and point mass with surface models with accretion/surface 
radii of 0.1 /?h- Throughout the low mass regime both types of 
model lead to results in good agreement with the analytic model 
of Type I migration. At and beyond a Saturn mass, where Type II 
migration becomes effective, the point mass with surface calcula- 
tions show slightly faster factor of 2) rates of migration. This 



is a result of the less evacuated gaps formed in sink with surface 
models, which results in a faster migration mechanism which is 
intermediate between Type I and II. The gap is more thoroughly 
evacuated when using an Accreting sink because it establishes no 
envelope structure to hinder accretion, and a negative pressure gra- 
dient to encourage it. This strips the corotational region of material 
quickly, resulting in a clearer gap at an equivalent time. However, 
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in both cases the gap is becoming more evacuated with time, which 
would be expected to slow migration to a Type II rate; this is a 
result of the models relatively short evolution of just 50 orbital pe- 
riods starting from a uniform disc. 

3.5 Protoplanets with surfaces and radiation hydrodynamics 

The structure of protoplanetary discs, and circumplanetary en- 
velopes and discs, are the result of competing gravitational and 
thermal effects. In the locally-isothermal models discussed to this 
point, the thermal effects associated with embedding a protoplanet 
in a protoplanetary disc have not been accounted for. To address 
this shortcoming we finally introduce radiative transfer and a real- 
istic equation of state into the migration models. The inclusion of 
a better thermodynamic treatment has several potentially important 
ramifications. The accretion energy released near a growing pro- 
toplanet warms its surroundings, which in concert with disc self- 
gravity leads to the development of self-consistent local disc and 
envelope structures. For example, a more poorly defined disc gap 
(opened by high mass planets) develops when radiative transfer is 
included (Fouchet & Mayer 2008j [AyUffe & Bate 2009a). Com- 
pressional heating also tends to smear out the spiral density waves 
that are launched by planets embedded in the discs, giving them 
lower peak densities than in locally-isothermal models. 

The impact of these changes can be seen in the panels of 
Fig. [To] which illustrate the torques acting on a 333 M© protoplanet 
in both a locally-isothermal model (left panel) and an interstel- 
lar opacity radiation hydrodynamical (RHD) model (right panel). 
The peak torques acting near the protoplanet are lower in the RHD 
case, though the maxima occur in similar locations (near the Roche 
lobe), illustrating the impact of the reduced local density. In the 
locally-isothermal case, the azimuthal angle at which torques from 
the corotation region fall away is considerably smaller than in the 
RHD case. This is a result of the more poorly cleared gap when us- 
ing RHD. Finally, the thermal spreading of the spiral arms can be 
seen by the broadened regions from which the associated torques 
are exerted, though there is no discernible difference in the mag- 
nitude of the torque acting from wave crest. Densities along a line 
in the midplane, 10 degrees clockwise of the same 333 M© proto- 
planet in both locally-isothermal, and RHD models, are shown in 
the right hand panel of Fig. [13] This distribution directly illustrates 
the difference in the vacuity of the gap. It is also clear that the den- 
sity maximum of the (outer) spiral arm is similar in both instances, 
despite the breadth being altered by the heating. The left panel is 
an equivalent plot for a 100 M© protoplanet. In this case the spiral 
arm density maximum is reduced by about a factor of 2 compared 
with the locally-isothermal case, possibly weakening the differen- 
tial Lindblad torques. However, the disc gap contains more mass 
than in the locally-isothermal model, once again leading to a mi- 
gration process even more similar to Type I, and so faster, than the 
expected Type II. The overall impact of using radiation hydrody- 
namics in this case will depend on the relative impact of these two 
opposing effects. 

The migration rates found in these interstellar grain opacity ra- 
diation hydrodynamical models are compared with those from oth- 
erwise equivalent locally-isothermal models in Fig. ^4] In the high 
mass regime, the impact is very small, with migration rates altered 
by no more than 30%. The 333 M© protoplanet, which maintains 
the same spiral arm densities but a less evacuated gap, migrates 
faster with RHD. This is consistent with an intermediate migration 
process, between Types I and II, as a result of the less evacuated gap 
formed with RHD. Once again, it is important to note that for both 
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Figure 12. The migration timescales of protoplanets modelled using Ac- 
creting sinks (asterisks) compared with those obtained using point masses 
with surfaces (plus signs), with accretion/surface radii of 0.1 Ru. These 
timescales are from models of locally-isothermal self-gravitating discs. The 
analytic Type I model of Tanaka et al. (2002) is plotted as a dashed line. 
With small accretion radii, the Accreting sinks give results inline with the 
analytic expectation for Type I migration and with the more realistic point 
mass with surface models. Beyond a Saturn mass, in the Type II regime, the 
Accreting sinks show slower migration by factors of 2. 



high mass protoplanets the change in migration timescales when 
radiative transfer is introduced are very small. In the lower mass 
Type I regime, the introduction of radiative transfer leads to sub- 
stantial increases in the migration timescales. A 50 M©protoplanet 
takes almost 4 times longer to migrate into its central star when 
modelled with radiative transfer than it does without, whilst a 33 
M© protoplanet takes more than 10 times longer. A 10 M© pro- 
toplanet is slowed by more than a factor of 2. In these cases, the 
densities around the protoplanets out to several /?h are lower when 
modelled with RHD using interstellar grain opacities than in the 
locally-isothermal models; about a 33 M© protoplanet the gas den- 
sity at the disc midplane at a radius of /?h is a factor of 3 lower. Sim- 
ilar reductions in Type I migration rates have been found in other 
recent work considering non-isothermal discs (e.g. | Paardekooper| 
& Mellema 2006). Recent models suggest that the slowing results 
from the dependence of the co-orbital torques on the radial entropy 
gradient of the protoplanetary disc, and that under some circum- 
stances these torques can be sufficient to reverse the direction of 
migration ( Paardekooper & Papaloizou|2 008a[|Bamteau & Masset 



2008b, Kley & Crida 2008, Kley et al. 2009, Paardekooper et al. 



|2010^ ). The reduced migration rates for low mass protoplanets are 
an important result in trying to increase the survival probability of 
proto-giant cores. Specifically it may justify the reduction, or some 
fraction of the reduction in Type I migration rates assumed in core 
accretion models that seek to ensure planets survive ( [Alibert et al.| 
|2004l|2005l[lda & Lin|2008l|Mordasini et al.|2009| ). 

In the point mass with surface models the envelope/disc sur- 
rounding the protoplanets is structured realistically due to the influ- 
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Figure 13. Midplane densities along lines 10 degrees clockwise of 100 M© (left) and 333 M© (right) protoplanets, modelled with point masses with surfaces of 
0.03 Ry{ after 50 orbits. The profiles run over rp ± 0.25rp. The solid lines mark the profiles for locally-isothermal calculations, whilst the dashed lines are taken 
from radiation hydrodynamics calculations using interstellar grain opacities. In both cases, the introduction of radiative transfer leads to poorer evacuation of 
the corotation region, allowing for greater torques to act from these radii. For the 100 M© protoplanet the outer spiral arm can be seen to have a lower peak 
density, a result of the heat produced in the shock. 
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Figure 14. The migration timescales of protoplanets modelled using point 
masses with surfaces of radius 0.03 in locally-isothermal (asterisks) and 
interstellar grain opacity radiation hydrodynamics (plus symbols) models. 
The analytic Type I model of Tanaka et al. (2002) is plotted as a dashed 
line. In the mass regime corresponding to Type I migration, the inclusion 
of radiative transfer with high opacities leads to slower rates of migration. 
For a 33 M© protoplanet the migration timescale is increased by an order of 
magnitude. Such slowed migration may improve the chances of proto-giant 
cores surviving long enough to grow to giant planet masses. 
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Figure 15. Scaleheights measured in discs containing a 33 M© protoplanet 
modelled using a point mass with surface of radius 0.03 Ru. The results 
of a locally-isothermal model are shown using a solid line, along with the 
radiation hydrodynamics models using 100 per cent (long-dashed), 10 per 
cent (short-dashed), and 1 per cent (dotted) interstellar grain opacities. The 
locally-isothermal and reduced opacity radiation hydrodynamics models all 
yield similar scaleheights through the disc, but the full opacity model gives 
a thicker disc, most notably around the radius of the embedded protoplanet. 
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Figure 16. Migration timescales for protoplanets in radiation hydrodynam- 
ics calculations, modelled by point masses with surface radii of 0.03 Ru, 
with three different opacities: 100 per cent (plus symbols), 10 per cent 
(asterisks), and 1 per cent (diamonds) interstellar grain opacities. The re- 
duced opacity models yield shorter migration timescales than their locally- 
isothermal equivalent models (c.f. Fig.[T4), with only the full opacity case 
leading to slower migration. Increasing the grain opacity of the disc leads 
to slower migration rates, with the effect becoming smaller for higher mass 
protoplanets. The analytic Type I model of |Tanaka et al., ( ,2002j is plotted as 
a dashed line. 



peaked. Fig. [TT] shows midplane densities and velocity vectors in 
the vicinity of the protoplanets. The velocity vectors make plain the 
spiral shocks, horseshoe orbits, and the transition to circumplane- 
tary flow deep within the Hill sphere. [20] shows density weighted 
temperature maps for the same models but concentrating on just the 
area surrounding the protoplanet. 

The diminution of the propagating density waves with increas- 
ing opacity maybe in part due to the process of wave channelling. In 
a vertically isothermal disc, these density perturbations propagate 
as plane waves (Lubow & Pringle_1993j . However, once radiative 
transfer is introduced, a vertical temperature profile is established, 
in which the midplane of the disc is hottest due to the effects of 
compression and viscous heating. This leads to wave channelling, 
causing the energy of the wave to quickly become confined to the 



surfaces of the disc, where it is dissipated ( Korycansky & Pringle 



ences of gravity, self-gravity, and thermodynamics. It is therefore 
possible to account for the differentiation of the angular momen- 
tum of accreted material into the planet's orbital angular momen- 
tum, and into rotation about the planet. This differentiation has next 
to no impact upon the compound protoplanet' s orbital angular mo- 
mentum. Typically a fraction of ~ 10~^ of the angular momentum 
of the gas accreted into a circumplanetary disc is diverted into the 
rotation of that disc about the planet. 

The impact of radiative transfer may go beyond the immedi- 
ate vicinity of a protoplanet, for example, inflating a protoplanetary 
disc's scaleheight as a result of viscous heating at the midplane. 
Thicker discs result in smaller torques from the inner and outer 
Lindblad resonances, scaling as (H/R)'^ ( Ward 1997 ), whilst the 
mismatch between the torques scales as H/R. This means that the 
migration rates scale as (H/R)'^, and a hotter, thicker disc should 
lead to slower migration ('Mas set & Kley|2066p . Fig. p3] shows the 
measured value of H/R against radius for a disc with an embedded 
33 M© protoplanet. It can be seen that there is a larger scaleheight 
at and about the protoplanet' s orbital radius in a radiation hydro- 
dynamics calculation using interstellar opacities when compared 
with all the other calculations. The scaleheight is a factor of ^ 1.2 
greater, which can only account for an increase in the migration 
timescale of about 50 per cent. 

3.5.1 The impact of opacity in RHD models 

Our final results come from varying the protoplanetary disc's opac- 
ity to determine its effect upon migration. Fig. [16] shows that the 
effect of varying the opacity on the migration rates is very small 
for the high mass planets. For the low mass planets there is a trend 
of increasing migration timescale with increased opacity. Interest- 
ingly, protoplanets in the low opacity RHD models migrate faster 
than in the locally-isothermal models. 

The largest change in the migration timescales of low mass 
protoplanets is seen when employing the interstellar grain opaci- 
ties, with the change between the 1% and 10% models, and the 
locally-isothermal models being somewhat smaller. In all but the 
interstellar grain opacity calculations, the densities around a given 
protoplanet are almost identical, and so the local torques should 
be similar. The spiral arms are also likely to be most broadened, 
and their peak density most reduced, in the highest opacity cal- 
culations, with lower opacities leading to structures more similar 
to those found in the locally-isothermal models. These structures, 
and maximum local torques can be seen in the torque maps shown 
in Fig. [it] which are produced for a 33 M© protoplanet modelled 
using a point mass with a surface radius of 0.03 Ru- The torques 
acting at and around the Roche lobe can be seen to be weaker in 
the interstellar opacity case (right most panel), than in any of the 
other cases, and the spiral arms most diminished. 

The diminution of the spiral arms with increasing opacity can 
be further seen in Figs. [18] and [19] which are surface density and 
density weighted temperature maps respectively of discs with em- 
bedded protoplanets, once again modelled by point masses with 
surfaces of radius 0.03 Ru- From left to right the protoplanet masses 
are 33, 100, and 333 M© respectively, and from top to bottom the 
models are locally-isothermal, and RHD with 1, 10, and 100 per 
cent interstellar grain opacities respectively; we have omitted the 
10 M© cases because the disturbance caused by the protoplanet is 
barely visible. Moving from the locally-isothermal density distribu- 
tion (top-row of figure [18]) to the interstellar opacity case (bottom 
row) it is evident that where a gap forms, it is less clear at higher 
opacities, and that the density waves propagating outwards are less 



1995 


ILubow & Ogilvie|1998||Ogilvie & Lubow|1999||Bate et al. 


2002 


1. This results in the more limited propagation of the spiral 



arms seen in the non-isothermal models, and which is clearest for 
low mass planets where the weak spirals barely last an entire plan- 
etary orbital period. It is also possible to see in Fig. [TT] that the 
breadth of the spiral arms increases with opacity, as was discussed 
earlier. 



In Fig. [19] and Fig. [20] it is possible to see the hot gas 
that develops around the protoplanets, particularly those of higher 
mass, when modelled with RHD and interstellar grain opacities. 
In Fig. [19] it is also possible to see the heat produced in the spiral 
shocks. The spatial extent of the hot gas surrounding a protoplanet 
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Figure 17. Torque distributions around a 33 M© protoplanet modelled by a point mass with a surface of radius 0.03 Ru after 50 orbits. From left to right the 
models are locally-isothermal, and then RHD with 1 per cent, 10 per cent, and 100 per cent interstellar grain opacity respectively, all including self-gravity. 
The Roche lobe is marked on in each case using a white line. All the torques above the white dashed line (y = 0) are positive, whilst all those below are 
negative. The spiral arms are particularly poorly defined in the interstellar grain opacity case (right most panel), and the torques acting from around the Roche 
lobe are weakest in this case. 



reduces with the reduction of the opacity, which allows the less 
dense outer reaches to cool more easily. It is also possible to see 
the higher temperatures in the spiral shocks, even in those waves 
that have propagated far from the planet, when compared with the 
underlying temperature distribution shown in the top row. 



4 DISCUSSION 

It seems that when planets form, they should also migrate. This 
migration poses problems for the survival of gas giant cores. Con- 
versely, such migration may aid in planet growth by opening up 
undepleted regions of the disc for these protoplanets to accrete ( Al- 
[ibert et aL]2004| . It is known from the ever increasing number of 
exoplanet detections that a large population of planets do grow and 
survive around many different types of star. The work discussed 
in this paper, like the work of many authors, was undertaken to 
try and ascertain the influence of different physical processes upon 
migration. The models we present in section [33?T] are the most self- 
consistent to date, including self-gravity, radiative transfer, and pro- 
toplanets modelled with surfaces to produce realistic envelopes and 
circumplanetary discs. 

4.1 Code comparison 

We first made comparisons with the results of a grid code to ensure 
that our SPH model was suited to tackling the problem of migra- 
tion, a capacity in which it has been seldom used before. The abil- 
ity of the SPH method to model protoplanet migration was consid- 
ered by |de Val-Borro et aL] ( |2006| ), who made comparisons between 



many grid codes and two SPH codes. They found reasonable agree- 
ment in disc structures around embedded protoplanets, though the 
SPH calculations under-resolved some disc features. This differ- 
ence is mitigated to an extent in our models by the increased num- 
ber of particles, more than 6 times greater than the number used 
in these previous comparisons. Our migration rates obtained us- 
ing non-self-gravitating locally-isothermal calculations are in good 
agreement with the analytic model of| Tanaka et al.|p002| ), and 
so also with many previous similar numerical models which show 
similar agreement to the analytic model. 

4.2 Self-gravity 

We introduced disc self-gravity to examine its impact upon mi- 
gration. |Pierens & Hure| { [2005 1 examined the inclusion of disc 
self-gravity analytically. They broke down the analysis into two 
competing influences, that of the circumstellar disc-planet inter- 
action, and the circumstellar disc-circumplanetary disc interaction 
(gas self-gravity). All of our SPH models include the first of these 
two interactions, which [Pierens & Hure| ( |2005 ) found to acceler- 
ate Type I migration. So by introducing self-gravity, we mean the 
addition of the second interaction. IPierens & Hure found this to 
asymmetrically shift the outer and inner Lindblad resonances such 
that the migration rate reduces. More recently, numerical models 
were performed by Baruteau & Masset] ( |2008b l) to further explore 
the impact of including self-gravity. They ran models to explore the 
impact of both interactions discussed here, and found, in agreement 
with^Pieren s & Hure ( 2005) and the previous numerical models of 
[Nelson & B enz (2003), that the disc-disc interaction reduced the 
Type I migration rates. The consistent, though small, increases in 
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migration timescales found for the SPH calculations including self- 
gravity (Fig.|7]) are in agreement with these previous results. Whilst 
the effect is small in our models, in discs with masses several times 
the minimum mass solar nebula self-gravity might slow migration 
by a more considerable factor ( [Baruteau & Masset|2008b| ). 

4.3 Radiation hydrodynamics 

Introducing radiative transfer has allowed us to investigate migra- 
tion beyond the more common locally-isothermal models. Previous 
studies concerned with the impact of non-isothermality, have found 
that the inclusion of more complicated thermodynamics leads to in- 
creased migration timescales ( |Jang-Condell & Sa sselov 20 05]), and 
even the reversal of the Type I migration dir ection (,Paardekooper| 
|& Mellema|2006l|2008b||Kley & Crida|2008l|Kley et al.|2009| ). We 
are also able to obtain outward migration for a 10 M© protoplanet 
when modelling it with an Accreting sink particle of radius 0.03 /?h 
with radiation hydrodynamics in an interstellar grain opacity disc, 
as shown in Fig. [21] 

However, once we progress to using our point mass with sur- 
face model, which is the most realistic case, we find no cases of 
outward migration. The 10 M© case modelled in this way, with a 
surface radius of 0.03 /?h and using interstellar grain opacities, is 
shown in Fig. [21] for comparison with the Accreting sink model. 
[Kley et ar] ([2009 ) find outward migration to result from a density 
peak in the planet's immediate vicinity, just ahead of it, which pulls 
it around in its orbit, increasing its angular momentum and so caus- 
ing it to migrate outwards. To test the reliability of this result they 
experiment with different resolutions and various gravitational soft- 
ening lengths, the shortest being 0.5 /?h, to explore their impact on 
the resulting migration. The resulting outward migration was found 
to be invariant to these changes. The migration rates are based on 
calculations of torques acting on planet at a fixed orbital radius of 
5.2 AU, where the torques are scaled depending upon their distance 
from the protoplanet, d, according to a fraction obtained from the 
tapering function 



fbid) : 



exp 



d/Rn - b 
b/lO 



+ 1 



(5) 



This function is such that at b, the torque cut-off radius in units of 
Ru, the torques are halved (where they typically use b = 0.8). This 
tapering is applied to avoid large and noisy torques from the poorly 
resolved Hill radius where the resolution is just 3.3 gridcells/i?H- 
Changing the value of b to 0.6 altered their total torques by 10 
per cent in locally-isothermal models, and 30 per cent in the radi- 
ation hydrodynamics models, but the resolution limits their ability 
to usefully reduce the value of b further. In Fig. [TOj where our un- 
derlying model resolution is ^ 20 times higher than that of |Kley| 
|et al.| ( |2009| ), it can be seen that the greatest torques are acting 
from within the Roche lobe. In fact the peak torques are acting 
from a radius of ^ 0.9 Ru in the radiation hydrodynamics case 
where gas is still passing through the Hill sphere, and able to carry 
away angular momentum. The ability of the SPH model to resolve 
these torques clearly, removes the need to use a tapering function of 
the kind above, which otherwise would have scaled these torques 
down to 73 per cent of their true values. Furthermore, using our 
planet surface treatment yields a self-consistent envelope around a 
20 M© protoplanet with densities that are an order of magnitude 
higher than those obtained using gravitational softening in |Kley| 
|et al.[ s otherwise similar disc model. Such a change in density is 
significant given the source of the positive torques suggested by 
[Kley et aL] Given the improved resolution and realism of the planet 



treatment employed in our SPH models, our results suggest that 
outward migration is not an inevitable consequence of introducing 
more complex thermodynamics. 

We find that introducing radiative transfer to our most real- 
istic models of low mass protoplanets in high opacity discs can 
slow migration rates by up to an order of magnitude. The reduc- 
tion of Type I migration rates may be a result of enhanced corota- 
tion torques due to entropy gradients across the co-orbital region 
(fPaardekooper & Papaloizou'"2008a", Baruteau & Masset""2008b[ 
Kley & Crida 2008, Kley et al. 2009, Paardekooper et al. 20To|)! 
For high mass protoplanets, those capable of forming circumplan- 
etary discs { [Ayliffe & Bate|2009b| ), the dominance of gravity over 
thermal effects leaves the migration rates very little altered upon 
the introduction of radiative transfer when compared with locally- 
isothermal conditions. This competition between thermal and grav- 
itational influences was also found to play an important role in de- 
termining accretion rates onto growing protoplanets by [Ayliffe &] 
[Bate [ ( [2009a] ). As with migration, the accretion rates in RHD mod- 
els were found to be most dependent upon the grain opacity for 
low mass protoplanets, where the gravitational forces at work are 
smaller. [Ay Hffe & Bate|p009a] ) found that a 33 M© protoplanet, in 
an interstellar grain opacity disc, is expected to double in mass in 
^5x10"^ years. We find the migration timescale for the same pro- 
toplanet to be ^ 5x10^ years. Under the same conditions a 100 M© 
protoplanet will double in mass in ^ 9 X 10^ years, and has a mi- 
gration timescale of ^ 2 x 10"^ years. Whilst a 333 M© protoplanet 
doubles in mass in 4 X 10^ years and migrates into its central 
star in ^ 10"^ years. These results suggest that is should be common 
for a protoplanet embedded in an interstellar grain opacity disc to 
reach a Jupiter mass. However, to reach larger masses may require 
higher disc masses or some additional process to either slow a pro- 
toplanet' s migration, or accelerate its accretion. 

In the low mass regime the increase in migration timescales 
that we find using RHD may help remedy the problem of proto- 
giant core depletion that was identified by 'Alibert et al. (2004[), 
[Alibert et al. ( 2005] ), [Ida & L in (2008 ) and Mor dasini et al. (2009| . 
Ida & Lin] showed that the efficiency of Type I migration must be 
reduced by at least an order of magnitude to enable the formation 
of a population of planets that match current observed prevalences. 
We find that the inclusion of radiative transfer, in conjunction with 
a realistic treatment of accretion on to the protoplanet can provide 
just such a reduction in the low mass regime. 

As was discussed in the results section, in the Type II regime 
the introduction of radiative transfer has two significant effects 
which have opposing impacts on migration rates. The reduction 
in the peak density of the spiral arms due to heating results in 
lower torques from these features, slowing migration, whilst the 
poorer evacuation of the disc gap leads to a migration process that 
is intermediate between Types I and II. In all our high mass cases 
which are modelled with either Accreting sinks, or with surfaces, 
in locally-isothermal or radiation hydrodynamics models, the rate 
of migration is considerably faster than the discs viscous migration 
timescale. Only using the unrealistic Killing sinks do our Type II 
rates approach the analytic model for Type II migration of |Ward[ 
( [1997[ ), most likely because they can most rapidly clear a gap of 
material. Some component of the rapid high mass migration rates 
found here may result from the models starting with an unperturbed 
disc, which over 50 orbits fails to establish a gap of the depth that 
might be expected if the planet had evolved to its initial mass in 
situ. However, radiative transfer might well be expected to hinder 
the transition from Type I to the slower Type II migration due to in- 
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log column density [g/cnn^] 

Figure 18. Surface density maps of protoplanetary discs after 50 orbits of the embedded protoplanets. From left to right the protoplanet masses are 33, 100, 
and 333 M©. From top to bottom the models are locally-isothermal, and then radiation hydrodynamical with 1,10, and 100 per cent interstellar grain opacities. 



© 0000 RAS, MNRAS 000, 000-000 



1 8 B.A. Ayliffe & M.R. Bate 



1 ' 1 ' 1 

o 

1 ' 1 ' \ 


1 ' 1 ■ 1 




1 ■ 1 ■ 1 

o . 

1 ' 1 ' 1 


o 

1 1 1 1 1 


o 




1 1 1 1 1 


I I 




t 


I I 
I I 




10 AU 

I.I.I 


: » : 

I.I.I 


1 1 1 1 1 

I.I.I 





1.5 2 2.5 3 



log gas temperature [K] 

Figure 19. Density weighted temperature maps of protoplanetary discs after 50 orbits of the embedded protoplanets. From left to right the protoplanet masses 
are 33, 100, and 333 M®. From top to bottom the models are locally-isothermal, and then radiation hydrodynamical with 1, 10, and 100 per cent interstellar 
grain opacities. 
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Figure 20. Magnified sections of the whole disc midplane temperature maps shown in Fig. [19] focussing on the region surrounding the embedded protoplanets. 
Again, from left to right the protoplanet masses are 33, 100, and 333 M®. From top to bottom the models are locally-isothermal, and then radiation hydrody- 
namical with 1, 10, and 100 per cent interstellar grain opacities. It is particularly apparent for the 333 M© protoplanet that the energy released as a result of 
accretion leads to a larger region of hot gas, which is less able to cool as the opacity is increased. It is also possible to see the broadening of the spiral arms as 
the opacity is increased. 
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Figure 21. The orbital evolution of a 10 M© protoplanet modelled three 
ways, twice as an Accreting sink particle and once as a point mass with sur- 
face. The Accreting sink models are differentiated by the chosen equation of 
state, one locally-isothermal (marked isothermal above), and one using ra- 
diation hydrodynamics (marked RHD) which migrates outwards. The point 
mass with surface is modelled using RHD, and shows a slow but inwards 
migration. 

creased thermal pressure that acts to oppose gap formation, as seen 
in our models. 



locally-isothermal and RHD models. These accretion and migra- 
tion timescales suggest that a protoplanet embedded in an interstel- 
lar grain opacity disc might grow to a high mass within its available 
migration time. 

• We find that the introduction of radiative transfer to a model 
of a 10 M© protoplanet, represented by an Accreting sink particle 
embedded in an interstellar grain opacity disc, causes it to migrate 
outwards. However the inclusion of radiative transfer does not lead 
to outward migration when we employ a well resolved and more 
realistic model of the protoplanet with a surface. 
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5 SUMMARY AND CONCLUSIONS 

We have conducted models of protoplanet migration with a range 
of physics and a raft of protoplanet treatments. Beginning with 
locally-isothermal, non-self gravitating discs, we incrementally 
added physics to reach self-gravitating radiation hydrodynamic cal- 
culations. We progressed from protoplanets modelled by crude sink 
particles, to those modelled by gravitating point masses with plan- 
etary surfaces. The final models we produce are the most self- 
consistent to date, and turn up some interesting results. In partic- 
ular, 

• We establish that the 'inertial mass problem', suff'ered by pro- 
toplanets that develop circumplanetary discs in models without 
self-gravity, is an eff'ect of order 10 per cent for minimum mass 
models. 

• We find that Type I migration rates, those that aff'ect low mass 
planets (10-33 M©), can be slowed by up to an order of magni- 
tude by the inclusion of radiative transfer in protoplanetary discs 
with interstellar grain opacities. Such a reduction helps to justify 
the arbitrary reductions in the Type I migration rate used in planet 
synthesis models ( Alib ert et al.|2004| [20051 [Ida & Lin|2008lpor^ 
[dasini et"ar]|2009 ) where such a reduction is required to ensure 
planets survive beyond the disc lifetime. 

• The Type I migration timescale that we obtain for a 10 M© 
protoplanet, modelled using RHD in an interstellar grain opacity 
disc, is of similar order to its mass doubling time, as found in 
[Aylifi^e & Bate (2009a), whilst for a 33 M© protoplanet the migra- 
tion timescale is an order of magnitude longer than the mass dou- 
bling time. Similarly for a 100 M© protoplanet the mass doubling 
time is shorter than the migration timescale, but for such high mass 
protoplanets the migration timescale is largely unchanged between 
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